Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
Show the code
ggsave(sprintf("figs/%smain_figures/absolute_mortality.png", fig_suffix), width=8, height=5)relative_mortality_rate <- agg_df %>%mutate(mortality_rate = deaths / births *1000) %>%group_by(state, category) %>%mutate(mean_mr =mean(mortality_rate[time <"2022-03-01"])) %>%mutate(relative_mr = mortality_rate / mean_mr)relative_mortality_rate %>%filter(category =="total") %>%ggplot() +geom_smooth(aes(x=time, y=relative_mr, group=state, col=state), se=FALSE) +geom_point(aes(x=time, y=relative_mr, col=state), height=0, alpha=0.5) +theme_bw(base_size=16) +theme(legend.position =c(0.99, 0.99), legend.justification =c(1, 1),#legend.background = element_blank(), # Make legend background transparentlegend.title =element_blank() ) +scale_color_manual(values=c("red", "orange", "dark gray"),labels=c("Banned States (excl. Texas)", "Texas", "Non-banned States")) +ylab("Relative Infant Mortality Rate") +xlab("Year") +geom_vline(xintercept=lubridate::date("2022-01-01"), color="orange", linetype="dashed") +geom_vline(xintercept=lubridate::date("2023-01-01"), color="red", linetype="dashed") +scale_x_date(date_breaks ="1 year", # Set axis labels at yearly intervalsdate_labels ="%Y"# Format the labels to show the year )
Warning in geom_point(aes(x = time, y = relative_mr, col = state), height = 0,
: Ignoring unknown parameters: `height`
New names:
Rows: 1000 Columns: 15347
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," dbl
(15347): ...1, .chain, .draw, mu[1,1,1], mu[2,1,1], mu[3,1,1], mu[4,1,1]...
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`
[1] "total 2"
New names:
Rows: 1000 Columns: 3878
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," dbl
(3878): ...1, .chain, .draw, mu[1,1,1], mu[1,2,1], mu[1,3,1], mu[1,4,1],...
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`
[1] "congenital 3"
New names:
Rows: 1000 Columns: 7851
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," dbl
(7851): ...1, .chain, .draw, mu[1,1,1], mu[2,1,1], mu[1,2,1], mu[2,2,1],...
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`
[1] "neonatal 3"
New names:
Rows: 1000 Columns: 7851
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," dbl
(7851): ...1, .chain, .draw, mu[1,1,1], mu[2,1,1], mu[1,2,1], mu[2,2,1],...
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`
make_state_fit_plot(quantiles_df %>%filter(type=="total"), state_name="Ban States", category="Total", target="deaths") +theme_bw(base_size=16) +scale_x_date(date_breaks ="1 year", # Set axis labels at yearly intervalsdate_labels ="%Y"# Format the labels to show the year) +theme(plot.margin =margin(t =10, r =20, b =10, l =10) # Extend the right margin) +ggtitle("Model Fit - All states with bans")
Show the code
ggsave(sprintf("figs/%smain_figures/ban_states_fit_plot.png", fig_suffix), width=8, height=5)make_gap_plot(quantiles_df %>%filter(type=="total"), state_name="Ban States", category="Total", target="deaths") +theme_bw(base_size=16) +scale_x_date(date_breaks ="1 year", # Set axis labels at yearly intervalsdate_labels ="%Y"# Format the labels to show the year) +ggtitle("Gap Plot - All states with bans") +theme(plot.margin =margin(t =10, r =20, b =10, l =10) # Extend the right margin)
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `state = fct_relevel(state, "Ban States (excl. Texas)", "Ban
States")`.
Caused by warning:
! 1 unknown level in `f`: Ban States (excl. Texas)
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `state = fct_recode(state, `States w/ bans` = "Ban States",
`States w/ bans (excl. Texas)` = "Ban States (excl. Texas)")`.
Caused by warning:
! Unknown levels in `f`: Ban States (excl. Texas)
ftab <-make_mortality_table(all_samples, target_state="Ban States (excl. Texas)", tab_caption ="Expected difference† in infant deaths (count and rate) in all states excluding Texas that banned abortion in months affected by bans, overall, by race, timing of death, and cause of death.") # tab_footnote(html("Exposed months include January 2023 through December 2023 for the 13 states besides Texas that banned abortion at 6 weeks or completely.<br># † Expected differences are computed using the Bayesian hierarchical panel data model described in the methods section. We report expected counts and rates as the observed count (or rate) minus the expected difference.# The estimates for each row are computed using the same Bayesian model; note however that the expected percent change will not necessarily equal the percent change in expected(known as Jensen’s inequality).<br># * Categories for which the 95% credible intervals of the expected difference excludes zero."))ftab |>gtsave(sprintf("figs/%ssupplement_figures/ban_states_no_texas_table.png", fig_suffix), zoom=4)
# ftab <- make_mortality_table(all_samples, target_state="Texas")# ftab |> gtsave("~/Dropbox/abortion_results_and_data/mortality_figures/texas_table.png", zoom=4)state_and_treat_times <- all_samples %>%filter(exposure_code==1, state !="Ban States") %>%group_by(state) %>%summarize(treatment_start =min(time))treated_states <- state_and_treat_times %>%pull(state)treated_times <- state_and_treat_times %>%pull(treatment_start)for(i in1:length(treated_states)) { target_state <- treated_states[i]if(target_state =="Ban States (excl. Texas)"){ caption_string <-"all states with bans excluding Texas" }else{ caption_string <- target_state } tab_caption =sprintf("Supplementary Table %i. Expected difference† in infant deaths (count and rate) in %s in months affected by bans, overall, by race, timing of death, and cause of death.", i, caption_string)make_mortality_table(all_samples, target_state=target_state, tab_caption=tab_caption) |>gtsave(filename=sprintf("figs/%ssupplement_figures/tables/%s_table.png", fig_suffix, target_state), zoom=4)# tab_footnote(html("† Expected differences are computed using the Bayesian hierarchical panel data model described in the methods section. We report expected counts and rates as the observed count (or rate) minus the expected difference.# The estimates for each row are computed using the same Bayesian model; note however that the expected percent change will not necessarily equal the percent change in expected(known as Jensen’s inequality).<br># * Categories for which the 95% credible intervals of the expected difference excludes zero. <br># In the publicly available CDC Wonder data, cell counts between 1 and 9 are suppressed. Dash denotes these missing counts or rates; expected differences are still inferred via the Bayesian hierarchical model.")) |>make_state_fit_plot(quantiles_df %>%filter(type=="total"), state_name=target_state, category="Total", target="deaths") +theme_bw(base_size=16) +ggtitle(paste("Model Fit -", target_state, sep=" ")) +scale_x_date(date_breaks ="1 year", # Set axis labels at yearly intervalsdate_labels ="%Y"# Format the labels to show the year ) +theme(plot.margin =margin(t =10, r =20, b =10, l =10) # Extend the right margin )ggsave(filename=sprintf("figs/%ssupplement_figures/fit_and_gap_plots/%s_fit_plot.png", fig_suffix, target_state), width=8, height=5)make_gap_plot(quantiles_df %>%filter(type=="total"), state_name=target_state, category="Total", target="deaths") +theme_bw(base_size=16) +ggtitle(paste("Gap Plot -", target_state, sep=" ")) +scale_x_date(date_breaks ="1 year", # Set axis labels at yearly intervalsdate_labels ="%Y"# Format the labels to show the year ) +theme(plot.margin =margin(t =10, r =20, b =10, l =10) # Extend the right margin )ggsave(sprintf("figs/%ssupplement_figures/fit_and_gap_plots/%s_gap_plot.png", fig_suffix, target_state), width=8, height=5)}
black_states_tab <-make_state_table(all_samples %>%filter(type =="race"),target_category="Non-Hispanic Black",tab_caption ="Table 3. Expected difference† in non-hispanic black infant deaths (count and rate) split by states that banned abortion in months affected by bans.",footnote_text ="")# "Exposed months include January 2022 through December 2023 for Texas and January 2023 through December 2023 for other 13 states that banned abortion at 6 weeks or completely.<br># † Expected differences are computed using the Bayesian hierarchical panel data model described in the methods section. We report expected counts and rates as the observed count (or rate) minus the expected difference.# The estimates for each row are computed using the same Bayesian model; note however that the expected percent change will not necessarily equal the percent change in expected(known as Jensen’s inequality).<br># * Categories for which the 95% credible intervals of the expected difference excludes zero."black_states_tab |>gtsave(sprintf("figs/%smain_figures/states_table_black_deaths.png", fig_suffix), zoom=4)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: There were 7992 warnings in `summarise()`.
The first warning was:
ℹ In argument: `obs_sval = sqrt(...)`.
ℹ In group 1: `category = Hispanic` `.draw = 1`.
Caused by warning in `matrix()`:
! data length [735] is not a sub-multiple or multiple of the number of rows [22]
ℹ Run `dplyr::last_dplyr_warnings()` to see the 7991 remaining warnings.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: There were 1998 warnings in `summarise()`.
The first warning was:
ℹ In argument: `obs_sval = sqrt(...)`.
ℹ In group 1: `category = Total` `.draw = 1`.
Caused by warning in `matrix()`:
! data length [882] is not a sub-multiple or multiple of the number of rows [22]
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1997 remaining warnings.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: There were 5994 warnings in `summarise()`.
The first warning was:
ℹ In argument: `obs_sval = sqrt(...)`.
ℹ In group 1: `category = Congenital` `.draw = 1`.
Caused by warning in `matrix()`:
! data length [882] is not a sub-multiple or multiple of the number of rows [22]
ℹ Run `dplyr::last_dplyr_warnings()` to see the 5993 remaining warnings.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: There were 5994 warnings in `summarise()`.
The first warning was:
ℹ In argument: `obs_sval = sqrt(...)`.
ℹ In group 1: `category = Neonatal` `.draw = 1`.
Caused by warning in `matrix()`:
! data length [882] is not a sub-multiple or multiple of the number of rows [22]
ℹ Run `dplyr::last_dplyr_warnings()` to see the 5993 remaining warnings.
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: Removed 25000 rows containing missing values or values outside the scale range
(`geom_line()`).
Show the code
## Post-Hoc IMR for those who would have received abortion but did not## 478 excess deaths## 20,323 excess birthswrite_lines(sprintf("------ %i deaths over expectation divided by %i births over expectation * 1000 = %f Infant Mortality Rate for those seeking abortion who did not get one ------ ", 478, 20323, 478/20323*1000), "figs/main_figures/summaries_for_paper.txt", append=TRUE)### Variance Decompositioncov_dat <-read_csv("data/dobbscovariates_2024_02_07.csv")
Rows: 4 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): state
dbl (2): median_time_post_dobbs, median_time_pre_dobbs
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Call:
lm(formula = 100 * (causal_effect_ratio - 1) ~ category + driving_time_diff,
data = .)
Residuals:
ALL 4 residuals are 0: no residual degrees of freedom!
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.966 NaN NaN NaN
categoryNon-Hispanic Black 9.611 NaN NaN NaN
categoryNon-Hispanic White 6.653 NaN NaN NaN
categoryOther 5.753 NaN NaN NaN
driving_time_diff NA NA NA NA
Residual standard error: NaN on 0 degrees of freedom
(60 observations deleted due to missingness)
Multiple R-squared: 1, Adjusted R-squared: NaN
F-statistic: NaN on 3 and 0 DF, p-value: NA